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1 Introduction 

The practical use of realistic head models for 
solving the forward and inverse problems in 
EEG/MEG or ECG/MCG, is limited by the time- 
consuming amount of numerical calculations. Often, 
a tradeoff should thus be accepted between 
numerical accuracy and computing time. Here, we 
propose a general method allowing for the use of 
realistic models with much higher speed than 
currently available when many forward calculations 
have to be performed, and compatible with any kind 
of numerical methods (Boundary Element, Finite 
Element, Finite differences). It is here illustrated in 
the case of EEG and MEG, but can be applied to 
ECG and MCG as well. 

2 Methods 

For a given position r in the brain, the three Lead- 
Fields L Xi0 (r), L yJ) (r), and L z0 (r), corresponding to 
unit dipoles at location r with orientation along 
positive x, y, and z 3D space cartesian axis, 
respectively, form a basis of the space spanned by 
the forward solutions produced by all possible 
dipoles located at this location. Hence, any dipole at 
r with moment vector 

'fix' 

Q= Qy 

.fix. 

produces L(r) = {V(r),B(r)j satisfying : 

kf) = Qx L x,o( r )+QyL y .o(r) +Q z L z fi(r ). ( 1 ) 

Here, V(r) and B(r) stand for the vectors containing 
the potential and field values at all electrodes and 
sensors, produced by the dipole Q at r. 

2.1 Classical approach 

When using realistic models and the BEM, 
computing V(r) and B(r) amounts to solve the 
following linear systems 

v=K'c e v, ( 2 ) 

( 3 ) 


where : 

• V is the vector of potential values at each of the 
N nodes of the mesh, 

• A eeg is a N*N matrix, 

• V 0 is a vector depending on the source and the 
model geometry, 

• B is the vector of field components along the 
normal of each of the N c sensor coils, 

• B 0 is the vector of field components produced by 
the primary current alone which can be 
determined analytically, and 

• A meg is a N C *N matrix (A meg F being the 
contribution of volume currents). 

The matrices A eeg and A meg have to be computed 
only once, for they depend only on the geometry 
and conductivities of the different head 
compartments, and not on the sources. 

The resolution of the first linear system (Eq. 2) for 
the potential vector V is the most time-consumming 
part. Indeed, although the matrix A eeg is usually 
inverted or LU-decomposed as a preliminary step 
once for all, each computation of Eq. 2 still requires 
N 2 multiplications and N 2 additions, which makes 
the computation time still important and directly 
related to the mesh size. 

2.2 Proposed approach 

2.2.1 Preliminary computations 

Both the classical and the proposed approaches 
require the preliminary computation of the system 
matrices A eeg and A meg and the decomposition or 
inversion of A eeg . Our approach requires a 
subsequent preliminary computation, which has to 
be performed only once for a given subject and a 
given sensor configuration (electrodes and/or pickup 
coils). It should be noted that an experiment is often 
run once per subject with only one sensor 
configuration, so that this step should usually be 
done once per subject. 

This additional step consists in 1) the construction 
of a discrete grid covering the brain volume, and 2) 
the computation of the Lead-Fields L x 0 (r), L v „(r), 
and L z 0 (r), for each location r of the grid. 


B = B Q + A me9 V 




2.2.2 Estimation of the forward problem for an 
arbitrary source 

Once these preliminary step are completed, each 
forward problem is estimated from the precomputed 
Lead-Fields. Given an arbitrary dipole source A in 
the brain volume at position R d with moment vector 
Q d , the solutions Z, a = {V d , B d \ produced by this 
source are approximated as follows : first, the three 
Lead-Fields Z, Xi0 (i? a ), L yfi (R 3 ), and L z fR d ) are 
interpolated from the Lead-Fields corresponding to 
the grid sources neighboring A; second, L a 
corresponding to the arbitrary orientation of A is 
obtained through Eq. 1. Note that if A is a shallow 
source falling outside the grid, extrapolation is used. 
Three different interpolation techniques were 
considered here : the trilinear interpolation, the 
Bezier interpolation (Bernstein polynomials of order 
2), and the 3D spline interpolation [1,2], For the 3D 
spline interpolation, we used an order of 3 and the 
80 grid points nearest to the source. 

2,3 Evaluation procedure 

In these simulations, spherical geometry models 
were considered, and spherical isotropic grids were 
used with 6 different step sizes of 6, 8, 10, 12, 15, 
and 20 mm. The respective number of grid points 
were 9315, 3887, 2007, 1189, 619, and 251. 



Figure 1 : Electrode montage and sensor coil array 
used in the present simulations. 

First, to account for the intrinsic errors inherent to 
the interpolation approximation, the method was 
tested with analytically pre-computed Lead-Fields. 
Next, it was evaluated with Lead-Fields pre¬ 
computed using the linear BEM [3], Figure 1 shows 
the 1500-nodes-per-layer model, the 63-electrode 
EEG montage and the 143-channel MEG array 


(first-order gradiometer arrangement from CTF) 
used in these simulations. 

The method was evaluated for both forward and 
inverse problems, and both in EEG and in MEG. 

The forward problem accuracy was evaluated at 
each point of a 4.3-mm-step grid (25341 points), 
where we considered 3 unit dipole orientations in 
EEG (along positive x, y, and z axis), and 2 in MEG 
along the tangential vectors of the local spherical 
coordinate system (u e and u f ). Forward solutions 
interpolated from pre-computed Lead-Fields were 
compared to analytical solutions, using the error 
criterion : 
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where / Cj0r i stands for either the potential or the 
magnetic field value at the cth channel for a given 
orientation. Figure 2 gives an example of a 
precomputed Lead-Field grid used to interpolated 
Lead-Fields at each point of the fine grid. 



Figure 2 : In the forward problem evaluation, the 
Lead-Fields at each point of the fine grid (pluses) 
are interpolated from the Lead-Fields of a 10-mm- 
step grid (filled circles). 

Localization errors of the inverse problem were 
calculated for 270 dipole locations distributed along 
9 axis of the inner sphere upper right quadrant 
(figure 3), with 3 orientations along x, y, and z in 
the EEG case (total of 810 dipoles), and for two 


















tangential orientations in the MEG case (540 
dipoles). Because the meshed model was spherical, 
the dipole orientation was constrained to remain 
tangential during the inverse procedure in the MEG 
case. Localization and orientation errors were 
averaged across the 9 axis and the 2 or 3 
orientations to get mean errors for each depth. 



Figure 3: Position of the sources in the inverse 
evaluation. 

2.4 Effect of grid cropping 

Shallow sources produce high numerical errors on 
the forward problem, especially in EEG when no 
mesh refinement is used [4,5], However, the use of 
locally refined models is not very realistic in 
practice, since models should be modified as the 
source position changes during the inverse 
procedure, each time leading to the re-computation 
and inversion or decomposition of the system 
matrices A eeg and A meg . 

Using our proposed approach, shallow grid points 
should thus intuitively not be used. We tested 
whether cropped grids would bring any 
improvement on the accuracy of the interpolated 
forward solutions and the inverse localization 
procedure, by removing shallow grid points having 
a depth smaller than the mean length of all triangles 
(8 mm in our case). 

3 Results 

3.1 Computation time 

Preliminary computations associated to the 
construction of matrices A eeg and A meg , and to the 
decomposition or inversion of the matrix A eeg , are 
common to the classical and to the proposed 


approaches. For the simulations presented here, 
these computations took 61 minutes on a 333-Mhz 
Pentium II with 256 Mo RAM (using the LU- 
decomposition of the matrix A eeg ). The proposed 
approach lead to an additional preliminary 
computation corresponding to the pre-computed 
Lead-Fields over the grid, which was proportionnal 
to the number of grid points : about lh 7' in EEG 
and lh 17' in MEG for 1000 grid points. Figure 4 
shows the dramatic decrease in computation time 
when estimating all Lead-Fields of the fine grid 
using the proposed approach with respect to using 
the classical approach. 



Figure 4 : Computation time with all approaches. 


The trilinear interpolation lead to the highest gain 
factor in computation time of about 20000 in EEG 
and 13700 in MEG. The Bezier interpolation lead to 
gain factors of 1200 in EEG and 750 in MEG, and 
the 3D spline only to gain factors of 40 in EEG and 
23 in MEG. 

3.2 Intrinsic precision of the method 

Intrinsic forward RMS errors below 2% and inverse 
localization and orientation errors below 0.2 mm 
and 0.2 degrees could be easily achieved with a 
sufficiently small grid step (6-8 mm) and either 
Bezier or 3D spline interpolations for any dipole 
depth (figure 5). 

When deep dipole were considered (although not 
too near the center in the MEG case), larger grid 
steps (10-12 mm) could be used to obtain similar 
accuracy. The trilinear interpolation method, 
although faster than the 2 others, lead to less 
accurate results, especially in MEG and especially 
at the center of the model. 
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Figure 5: Intrinsic localization errors (mm) for 2 
different step sizes (6 & 8 mm) and the 3 different 
interpolation techniques. 

3.3 Comparison with the classical approach 

In case of BEM pre-computed Lead-Fields, grid 
steps of about 6-8 mm lead to very comparable 
localization errors than the classical BEM, 
irrespectively of the dipole depth. In other words, 
errors due to the interpolation are negligible with 
respect to numerical errors inherent to the BEM. For 
deep dipoles (depth > 20 mm), larger grid steps (10- 
12 mm) might even be used to obtain a good 
accuracy. In general, the Bezier interpolation 
showed superiority with respect to the trilinear 
interpolation and lead to more stable results than the 
3D spline interpolation. It should thus be prefered to 
the two others. 

3.4 Effect of grid cropping 

Grid cropping strongly decreased the orientation 
errors in the case of EEG for shallow sources 
(figure 6). This held for any interpolation method 
used and any grid step between 6 and 15 mm. Such 
improvement was only very slight in MEG, for 
which orientation errors were already very good for 
shallow sources. Localization errors were not 
improved by grid-cropping. 
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Figure 6 : Effect of grid cropping on the orientation 
errors in the EEG case. The grid step is 8 mm. Stars 
stands for classical BEM, squares circles and dots 
for for trilinear, Bezier, and 3D spline 
interpolations, respectively. 

4. Conclusion 

This method allows for the use of realistic models 
with speed comparable to those provided by the 
spherical model, without noticeable loss of accuracy 
(~ 0.2 mm) comparing to the precision of current 
EEG/MEG inverse methods (>lmm). It even 
increase of accuracy on dipole orientation 
estimation for shallow sources in the EEG case. 
Another advantage stems from the fact that only the 
lead-field precalculation time depends on the 
complexity of the model used, but not the 
interpolation procedure. 

A grid step of 8 mm and the Bezier interpolation 
were found to be good compromise between 
accuracy and Lead-Field precomputation time. 

It should be noted that the gain in computation time 
is inversely related to the number of sensors. 
However, as shown here with realistic sensor 
configurations, it brings appealing improvements 
even for currently used large sensor arrays. It would 
be especially useful in stochastic minimisation 
algorithm such as simulated annealing or genetic 
approaches, which require numerous forward 
calculations. 
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